% rc9 Feb-2011 - 
%  data quality checks - b_norm, e_zero, esum, e_avg & e_03d
%
% $Id: e_sum.m 488 2005-12-15 20:15:45Z edmaclin $
% e_sum.m - function version of eros_avg2.m; sum eros data per .sto timing info
%
% Ed Maclin eros_avg 12/11/02; e_avg 3/5/03
% $Log: e_sum.m,v $
% Revision 1.20  2004/09/07 19:10:26  odellj
% Cleaned up output.
%
% Revision 1.19  2004/09/07 18:35:33  odellj
% Fixed a minor conflict.
%
% Revision 1.18  2004/08/26 19:06:16  odellj
% Fixed with inputed out_path, in_path variables to allow for averaging of norm data.
%
% Revision 1.17  2004/08/25 22:12:03  odellj
% Added the ID tag and the ability to pass that into the id_tags structure if this routine is called.
%
% Revision 1.16  2004/08/18 18:32:24  edmaclin
% solved (?) failure to avg "norm" data (use "temp" filename)
%
% Revision 1.15  2004/08/16 19:13:59  odellj
% Code was halting if the fid didn't exist.  Created a conditional to account for that.
%
% Revision 1.14  2004/07/19 21:47:28  odellj
% Commiting from the previous, correct version.  Exist is not working under compilation.  (See 1.7 or the exist function branch.)
%
% Revision 1.13  2004/07/19 18:25:25  odellj
% if (exist([in_path file_name]) & exist(evt_file_name))
%
% --missing a '('
%
% Revision 1.12  2004/07/15 19:25:26  edmaclin
% cleaned up file checking code
%
% Revision 1.11  2004/07/15 17:46:56  edmaclin
% *** empty log message ***
%
% Revision 1.10  2004/07/15 17:41:43  edmaclin
% fclose 'all', not specific fids (which may not exist)
%
% Revision 1.9  2004/07/15 17:33:17  edmaclin
% moved 1st log print
%
% Revision 1.8  2004/07/15 15:38:38  odellj
% More error checking.
%
% Revision 1.7  2004/07/09 21:50:14  odellj
% Updated the logging of the subroutines.  Added log_file as a parameter for each of their calls.
%
% Revision 1.6  2004/07/07 18:07:53  odellj
% Prepped for logging.
%
% Revision 1.5  2004/06/30 18:13:48  odellj
% Added more verbose fprintf statements.
%
% Revision 1.4  2004/05/10 21:44:49  edmaclin
% declared SE arrays to avoid missing (n=0) bins
% process only "good" channels
%
% Revision 1.3  2004/03/31 21:52:04  edmaclin
% runs 1000 times faster ;-)
%
% Revision 1.2  2004/03/22 20:15:37  odellj
% Added optionality of .nir files
%
% Revision 1.1.1.1  2004/02/18 20:09:49  odellj
% initial import into CVS
%
% Revision 1.2  2004/02/05 20:42:52  odellj
% Added Id and Log tags.
% 

function e_sum(in_path,out_path,log_file)

global hdr file_name i_block dc ac ph adc aac aph n2 ssdc ssac ssph
global nir_ph nir_ac nir_dc msg zap_flg id_tags
global good_dc_pts good_dc_chs good_ac_pts good_ac_chs
global good_ph_pts good_ph_chs good_chs event n_good

fprintf('e_sum\n');

% in_path=[hdr.data_path hdr.exp_name hdr.data_type hdr.hi '-' hdr.lo '\'];
% out_path=[hdr.data_path hdr.exp_name hdr.data_type hdr.hi '-' hdr.lo '\'];

e_sum_id = '$Id: e_sum.m 488 2005-12-15 20:15:45Z edmaclin $';
check = strcmp(e_sum_id,id_tags); % Make sure we don't already have the tag
if sum(check) == 0
    id_tags{size(id_tags,1)+1,1} = e_sum_id;  % Add it in then.
end

scrsz = get(0,'ScreenSize'); % for debug plotting

% Determine whether input and output files exist
%  !!! 'exist' function doesn't work in compiler :-(
%
in_file_flg=fopen([in_path file_name]);
evt_file_name=[hdr.data_path hdr.exp_name 'evt\' file_name];
evt_file_flg=fopen([evt_file_name]);
k=findstr('.',file_name);
out_file_name=[file_name(1:k(1)) 'avm'];
nir_file_name=[file_name(1:k(1)) 'nir'];
xls_file_name=[file_name(1:k(1)) 'xlsx'];

out_file_flg=fopen([out_path out_file_name]);

if in_file_flg > 1, fclose(in_file_flg); end
if evt_file_flg > 1, fclose(evt_file_flg); end
if out_file_flg > 1, fclose(out_file_flg); end

if ~strcmp(log_file,'')
    logid = fopen(log_file,'a');
else logid = 1; end

fprintf(logid,'Summing...\n');

if out_file_flg < 0
	if ((in_file_flg > 0) & (evt_file_flg > 0))       
        % Load 'event' matrix in .mat format 
        %
        msg=sprintf('Reading %s',evt_file_name);
        load(evt_file_name,'event','-mat');
        fprintf(logid,'Reading %s\n',evt_file_name);

        % load EROS data in .mat format
        %
        temp=[in_path file_name];
        load(temp,'dc','ac','ph','mrdc','mrac','-mat');
        fprintf(logid,'Reading %s\n',temp);

		                    
        % Sum entire block (for hemodynamic & spectral analyses)
        % NOTE: May remove or make optional.  Not usually used.
        disp('here1')
        fprintf(logid,'Summing blocks...\n');
        disp('here2')
        nir_ph=nir_ph+ph;
        nir_ac=nir_ac+ac;
        nir_dc=nir_dc+dc;
        
        n_trials=size(event,1);
        
        good_chs=zeros(1,hdr.n_chans);
        
        % check each channel for mrac, mrdc
        % done here (not b_norm) to allow use of liberal "norming"
        %
        disp('Setting good block chans')
        good_mrac_chs=find(mrac>hdr.qual.mrac_thr);
        good_chs=good_mrac_chs;
        hdr.n_good.mrac(good_mrac_chs)=hdr.n_good.mrac(good_mrac_chs)+1;

        good_mrdc_chs=find(mrac>hdr.qual.mrdc_thr);
        good_chs=intersect(good_chs,good_mrdc_chs);
        hdr.n_good.mrdc(good_mrdc_chs)=hdr.n_good.mrdc(good_mrdc_chs)+1;
        
%         fprintf('%3.0f good_mrdc_chs\n',good_mrdc_chs)
%         fprintf('%3.0f good_mrac_chs\n',good_mrac_chs)

        % Compute sums per condition and channel
        %
        for i_trial=1:n_trials
            Summing blocks
            i_cond=event(i_trial,2); % only include trials for current condition
            start=event(i_trial,1)-hdr.avg_baseline+1;
            stop=event(i_trial,1)+hdr.epoch_length;
            if (start<1) | (stop>length(dc))
                fprintf(logid,'i_trial=%4.0f,  start=%4.0f,  stop=%4.0f, length(dc)=%4.0f\n',...
                    i_trial,start,stop,length(dc))
            else
                % check each channel for dc range, ac min, delay range
                %   & modulation min for this epoch
                %
                disp('Setting good trial chans')
                dc_max=max(dc(start:stop,:));
                dc_min=min(dc(start:stop,:));
                dc_range=dc_max-dc_min;
                good_dc_chs=find(dc_range<hdr.qual.dc_range);
                good_chs=intersect(good_chs,good_dc_chs);
                hdr.n_good.dc(good_dc_chs)=hdr.n_good.dc(good_dc_chs)+1;
                
                % check each channel for ac min
                %
%                 disp('checking AC min')
%                 [good_ac_pts,good_ac_chs]=find(ac(start:stop,:)>hdr.qual.ac_min);
%                 good_chs=intersect(good_chs,good_ac_chs);
%                 hdr.n_good.ac(good_ac_chs)=hdr.n_good.ac(good_ac_chs)+1;

                % check each channel for delay range
                %
                disp('checking delay range')
                ph_max=max(ph(start:stop,:));
                ph_min=min(ph(start:stop,:));
                ph_range=ph_max-ph_min;
                good_ph_chs=find(ph_range<hdr.qual.delay_range);
                good_chs=intersect(good_chs,good_ph_chs);
                hdr.n_good.ph(good_ph_chs)=hdr.n_good.ph(good_ph_chs)+1;

                % check each channel for modulation min for this epoch
                % low modulation may mean AC clipping
                % 1st compute modulation for this epoch
                %
                tmrdc=repmat(mrdc,stop-start+1,1);
                tmrac=repmat(mrac,stop-start+1,1);
                epc_mod=(tmrac.*(ac(start:stop,:)+1))./(tmrdc.*(dc(start:stop,:)+1));
                %
                % use 'min' or 'mean' modulation during epoch ?
                %
%                 min_mod=min(epc_mod); % minimum mod on each channel
%                 good_mod_chs=find(min_mod>hdr.qual.mod_thr);
                mean_mod=mean(epc_mod); % mean mod on each channel
                good_mod_chs=find(mean_mod>hdr.qual.mod_thr);
                
                good_chs=intersect(good_chs,good_mod_chs);
                hdr.n_good.mod(good_mod_chs)=hdr.n_good.mod(good_mod_chs)+1;
                
                disp(['size(good_chs)= ' num2str(length(good_chs))])
                hdr.n_good.tot(good_chs)=hdr.n_good.tot(good_chs)+1;
                
                % increment accepted channels for this bin
                %
                n2(good_chs,i_cond)=n2(good_chs,i_cond)+1;
                
                % sum good channels for this trial
                %
                adc(:,good_chs,i_cond)=adc(:,good_chs,i_cond)+dc(start:stop,good_chs);
                ssdc(:,good_chs,i_cond)=ssdc(:,good_chs,i_cond)+dc(start:stop,good_chs).^2;
                aac(:,good_chs,i_cond)=aac(:,good_chs,i_cond)+ac(start:stop,good_chs);
                ssac(:,good_chs,i_cond)=ssac(:,good_chs,i_cond)+ac(start:stop,good_chs).^2;          
                aph(:,good_chs,i_cond)=aph(:,good_chs,i_cond)+ph(start:stop,good_chs);
                ssph(:,good_chs,i_cond)=ssph(:,good_chs,i_cond)+ph(start:stop,good_chs).^2;   
                
                % plot/print this epoch for debugging
                %
%                 close all
                %
                % plot dc range
                %
%                 f1=figure('Position',[.05*.9*scrsz(3) .05*.9*scrsz(4) .9*scrsz(3) .9*scrsz(4)]);
%                 set(gcf,'Color','w')
%                 for i_ch=1:64   % plot just 64 chans
%                     subplot(8,8,i_ch)
%                     hold on
%                     plot(dc(start:stop,i_ch),'r')
%                     xlims=xlim;
%                     ylims=ylim;
%                     text(1,ylims(2),['range=' num2str(dc_range(i_ch),'%4.0f')],'FontSize',8)
%                     plot([xlims(1) xlims(2)],[dc_max(i_ch) dc_max(i_ch)],'k')
%                     plot([xlims(1) xlims(2)],[dc_min(i_ch) dc_min(i_ch)],'k')
%                 end
%                 for i_ch=good_dc_chs   % plot just good_dc chans
%                     if i_ch<65          % and only up to ch 64
%                         subplot(8,8,i_ch)
%                         hold on
%                         plot(dc(start:stop,i_ch),'b')
%                     end
%                 end
                %
                % plot ph (delay) range
                %
%                 f2=figure('Position',[.05*.9*scrsz(3) .05*.9*scrsz(4) .9*scrsz(3) .9*scrsz(4)]);
%                 set(gcf,'Color','w')
%                 for i_ch=1:64   % plot just 64 chans
%                     subplot(8,8,i_ch)
%                     hold on
%                     plot(ph(start:stop,i_ch),'r')
%                     xlims=xlim;
%                     ylims=ylim;
%                     text(1,ylims(2),['range=' num2str(ph_range(i_ch),'%4.0f')],'FontSize',8)
%                     plot([xlims(1) xlims(2)],[ph_max(i_ch) ph_max(i_ch)],'k')
%                     plot([xlims(1) xlims(2)],[ph_min(i_ch) ph_min(i_ch)],'k')
%                 end
%                 for i_ch=good_ph_chs   % plot just good_dc chans
%                     if i_ch<65          % and only up to ch 64
%                         subplot(8,8,i_ch)
%                         hold on
%                         plot(ph(start:stop,i_ch),'b')
%                     end
%                 end
                
                % print mean (min?) modulation
                %
%                 disp(['hdr.qual.mod_thr=' num2str(hdr.qual.mod_thr)])
%                 for i_ch=1:64
%                     if epc_mod(i_ch) < hdr.qual.mod_thr
%                         disp([num2str(i_ch) ')' num2str(epc_mod(i_ch),'%4.2f') ' BAD'])
%                     else
%                         disp([num2str(i_ch) ')' num2str(epc_mod(i_ch),'%4.2f') ' GOOD'])
%                     end
%                 end
                
%                 pause % only for debugging...
            end
        end % i_trial
	
        fprintf('hdr.n_good.mrdc = ');fprintf('%4.0f',hdr.n_good.mrdc);fprintf('\n')
        fprintf('hdr.n_good.mrac = ');fprintf('%4.0f',hdr.n_good.mrac);fprintf('\n')
        fprintf('hdr.n_good.dc   = ');fprintf('%4.0f',hdr.n_good.dc);fprintf('\n')
%         fprintf('hdr.n_good.ac   = ');fprintf('%4.0f',hdr.n_good.ac);fprintf('\n')
        fprintf('hdr.n_good.ph   = ');fprintf('%4.0f',hdr.n_good.ph);fprintf('\n')
        fprintf('hdr.n_good.mod  = ');fprintf('%4.0f',hdr.n_good.mod);fprintf('\n')
        fprintf('hdr.n_good.tot  = ');fprintf('%4.0f',hdr.n_good.tot);fprintf('\n')
        
        msg=strcat(msg,sprintf('; n-total =      '));
        msg=strcat(msg,sprintf('%5.0f',n2));
	
        % collect all n_good vals in one matrix and write to .xlsx file
        %
        n_good=[1:hdr.n_chans;hdr.n_good.mrdc;hdr.n_good.mrac;hdr.n_good.dc;...
            hdr.n_good.ph;hdr.n_good.mod;hdr.n_good.tot];
        c_good={n_good};
        n_good_labels={'Ch #';'mrdc';'mrac';'dc';'ph';'mod';'tot'};
        warning off MATLAB:xlswrite:AddSheet
        xlswrite([out_path xls_file_name],n_good_labels,['Block' num2str(i_block,'%02.0f')],'A1')
        xlswrite([out_path xls_file_name],n_good,['Block' num2str(i_block,'%02.0f')],'B1:IW7')

        
	else % if various files exist or don't...
        if (in_file_flg < 0)
            msg=sprintf('%s does not exist.\n',[in_path file_name]);
            fprintf(logid,'%s does not exist.\n',[in_path file_name]);
        end
        if (evt_file_flg < 0)
            msg=sprintf('%s does not exist.\n',evt_file_name);
            fprintf(logid,'%s does not exist.\n',evt_file_name);
        end     
	end
	
else    % Out_file exists
    msg=sprintf('%s already exists!\n',[out_path out_file_name]);
    fprintf(logid,'%s already exists!\n',[out_path out_file_name]);
end

if logid ~= 1, fclose(logid); end